#' =============================================================================
#' FILE: 05_discussion.R
#' DESCRIPTION:
#'   Analyzes social distance and affective polarization using CFA and
#'   visualizations. Produces tables and figures for the discussion.
#'
#'   Produces table A7 and figure 08.
#'
#' PACKAGES REQUIRED: pacman, tidyverse, lavaan, ggpubr
#'
#' OUTPUTS:
#'   - 03_figures/fg_08.jpg
#' =============================================================================

# Packages and dataset ---------------------------------------------------------

# Required packages
if (!require("pacman")) install.packages("pacman")
pacman::p_load(tidyverse, 
               lavaan, 
               ggpubr)

# Data
df <- readRDS("04_outputs/clean_dataset.rds")

# Social Distance variables: table A7 ------------------------------------------

# Recode variables function
recode_st <- function(var) {
  var <- trimws(var)
  new_var <- NA
  new_var[var == "Definitivamente no haría esto"] <- 0
  new_var[var == "Probablemente no haría esto"] <- 1
  new_var[var == "Probablemente haría esto"] <- 2
  new_var[var == "Definitivamente haría esto"] <- 3
  return(new_var)
}

df$sd_1_1 <- recode_st(df$social_distance_1_1)
df$sd_1_2 <- recode_st(df$social_distance_1_2)
df$sd_2_1 <- recode_st(df$social_distance_2_1)
df$sd_2_2 <- recode_st(df$social_distance_2_2)
df$sd_3_1 <- recode_st(df$social_distance_3_1)
df$sd_3_2 <- recode_st(df$social_distance_3_2)
df$sd_4_1 <- recode_st(df$social_distance_4_1)
df$sd_4_2 <- recode_st(df$social_distance_4_2)

## 2-Factor Model
m_sd <- '
sd_fuji =~ sd_1_1 + sd_2_1 + sd_3_1 + sd_4_1
sd_anti =~ sd_1_2 + sd_2_2 + sd_3_2 + sd_4_2
sd_fuji ~~ sd_anti
'
# Fit
m_sd_fit <- sem(m_sd,
                missing = "FIML", # for missing values
                data = df)
## Summary and fit: Table A7
fitMeasures(m_sd_fit, c("cfi", "rmsea", "srmr")) 

df <- cbind(df, lavPredict(m_sd_fit))

# Social distance plot (Panel a) -----------------------------------------------
  
sd <- pivot_longer(df %>% 
                     select(cha_id, sd_fuji, sd_anti), 
                   cols = c(sd_fuji, sd_anti), names_to = "id", values_to = "sd") %>% 
  filter(cha_id != 'none', 
         !is.na(sd)) %>% 
  mutate(id = recode(id, 
                     "sd_anti" = "Social Distance Towards Anti-Fujimoristas",
                     "sd_fuji" = "Social Distance Towards Fujimoristas"), 
         id = factor(id, levels = c("Social Distance Towards Fujimoristas",
                                    "Social Distance Towards Anti-Fujimoristas"))) %>% 
  ggplot(aes(x = sd, color = cha_id)) +
  geom_density(alpha = 0.3) +
  scale_color_manual(values = c("grey", "black")) +
  scale_fill_manual(values = c("grey", "black")) +
  labs(y = "Density", x = "Value", color = "Identity", fill = "Identity") +
  theme(legend.position = "bottom", 
        legend.title.position = "top",
        legend.title = element_text(hjust = 0.5),
        legend.key = element_rect(fill="white"),
        legend.box.background = element_rect(colour = "black"),
        legend.background = element_blank(),  
        panel.background = element_rect(fill='white', colour='white'),
        panel.grid.major = element_line(colour = "grey95"),
        panel.grid.minor = element_line(colour = "grey95"),
        axis.ticks.y = element_blank(),
        axis.ticks.x = element_blank(),
        strip.background = element_rect(fill="white")) +
  facet_grid(~ id)
sd

# Feeling Thermometers (Panel b) -----------------------------------------------

# Ingroup and outgrup
df$ft_fuji <- as.numeric(df$pb02_9) - 50
df$ft_anti <- as.numeric(df$pb02_11) - 50

# Point estimates
point_est <- df %>% 
  select(cha_id, ft_fuji, ft_anti) %>% 
  pivot_longer(-cha_id) %>% 
  filter(!is.na(value) &
           cha_id != 'none') %>% 
  mutate(inout = case_when(
    cha_id == "Fujimorismo" & name == "ft_fuji" ~ 'Ingroup Favoritism',
    cha_id == "Fujimorismo" & name == "ft_anti" ~ 'Outgroup Prejudice',
    cha_id == "Anti-Fujimorismo" & name != "ft_fuji" ~ 'Ingroup Favoritism',
    cha_id == "Anti-Fujimorismo" & name != "ft_anti" ~ 'Outgroup Prejudice'),
    value = abs(value),
    cha_id = factor(cha_id, levels = c("Fujimorismo", 
                                       "Anti-Fujimorismo"))) %>% 
  group_by(cha_id, inout) %>% 
  summarise(scr_value = mean(value), 
            upr_value = mean(value)+1.96*(sd(value)/sqrt(length(value))),
            lwr_value = mean(value)-1.96*(sd(value)/sqrt(length(value)))) 

# Plot
ft <- df %>% 
  select(cha_id, ft_fuji, ft_anti) %>% 
  pivot_longer(-cha_id) %>% 
  filter(!is.na(value) &
           cha_id != 'none') %>% 
  mutate(inout = case_when(
    cha_id == "Fujimorismo" & name == "ft_fuji" ~ 'Ingroup Favoritism',
    cha_id == "Fujimorismo" & name == "ft_anti" ~ 'Outgroup Prejudice',
    cha_id == "Anti-Fujimorismo" & name != "ft_fuji" ~ 'Ingroup Favoritism',
    cha_id == "Anti-Fujimorismo" & name != "ft_anti" ~ 'Outgroup Prejudice'),
    value = abs(value), 
    cha_id = factor(cha_id, levels = c("Fujimorismo", 
                                       "Anti-Fujimorismo"))) %>% 
  ggplot() +
  geom_density(aes(x = value, color = inout), alpha = 0.15) +
  geom_point(point_est, 
             mapping = aes(x = scr_value, y = 0.04, color = inout), 
             position = position_dodge(width = 0.05)) +
  geom_errorbarh(point_est, 
                mapping = aes(xmin = lwr_value, xmax = upr_value, y = 0.04, color = inout), 
                position =position_dodge(width = 0.05),
                height = 0.006) +
  scale_color_manual(values = c("grey", "black")) +
  labs(x = "Score", y = "Density", color = "Affects") +
  theme(legend.position = "bottom", 
        legend.title.position = "top",
        legend.title = element_text(hjust = 0.5),
        legend.key=element_rect(fill="white"),
        legend.box.background = element_rect(colour = "black"),
        legend.background = element_blank(),
        panel.background = element_rect(fill='white', colour='white'),
        panel.grid.major = element_line(colour = "grey95"),
        panel.grid.minor = element_line(colour = "grey95"),
        axis.ticks.y = element_blank(),
        axis.ticks.x = element_blank(),
        strip.background = element_rect(fill="white")) + 
  facet_wrap(~ cha_id, 
             ncol = 2)
ft

# Combine Plot -----------------------------------------------------------------

ggarrange(
  sd,
  ft,
  nrow = 2, 
  labels = c("a)", "b)"),
  heights = c(1, 1)  # Adjust these values to change the relative heights of the rows
)

ggsave('03_figures/fg_08.jpg',
       height = 8,
       width = 9)